/*
 *  MatrixUtils.h
 *  NSE_SOLVER
 *
 *  Created by Patrick Zulian on 1/14/11.
 *  Copyright 2011 ICS - USI - Lugano. All rights reserved.
 *
 */

#ifndef NSE_SOLVER_MATRIX_UTILS_H
#define NSE_SOLVER_MATRIX_UTILS_H

#include <iostream>
#include <fstream>
#include <string>
#include <exception>
#include <math.h>
#include <vector>

#include "Base.h"
#include "Utils.h"

namespace nsesolver {
	
	template< class MatrixType >
	int NonZeros(const MatrixType &m)
	{
		int nnz = 0;
		for (int i = 0; i < m.size(); i++) {
			nnz += m[i] != 0;
		}
		return nnz;
	}
	
	template< class MatrixType >
	const MatrixType RowMean(const MatrixType &m)
	{
		MatrixType result(SumRows(m));
		result *= 1.0/m.columns(); 
		return result;
	}
	
	template< class MatrixType >
	const MatrixType ColumnMean(const MatrixType &m)
	{
		MatrixType result(SumColumns(m));
		result *= 1.0/m.rows(); 
		return result;
	}
	
	template< class MatrixType >
	const MatrixType AroundMean(const MatrixType &m)
	{
		MatrixType result(m);
		MatrixType mean(RowMean(m));

		for (int j = 0; j < result.columns(); j++) {
			for (int i = 0; i < result.rows(); i++) {
				result(i, j) -= mean(i, 0);
			}
		}
		return result;
	}
	
	template< class MatrixType >
	const MatrixType Diag(const MatrixType &m)
	{
		using std::min;
		const int n = min(m.rows(), m.columns()); 
		MatrixType result(n, 1);
		
		for (int i = 0; i < n; i++) {
				result(i, 0) = m(i, i);
		}
		return result;
	}
	
	
	template< class MatrixType >
	const MatrixType Cov(const MatrixType &m)
	{
		MatrixType temp(AroundMean(m));
		MatrixType result = (1.0/(temp.columns() -1.0)) * (temp * temp.transposed());
		return result;
	}
	
	
	template< class MatrixType >
	const MatrixType SumRows(const MatrixType &m)
	{
		MatrixType result(m.rows(), 1);
		for (int i = 0; i < m.rows(); i++) {
			result(i, 0) = 0;
			for (int j = 0; j < m.columns(); j++) {
				result(i, 0) += m(i, j);
			}
		}
		return result;
	}
	
	template< class MatrixType >
	const MatrixType SumColumns(const MatrixType &m)
	{
		MatrixType result(1, m.columns());
		for (int j = 0; j < m.columns(); j++) {
			result(0, j) = 0;
			for (int i = 0; i < m.rows(); i++) {
				result(0, j) += m(i, j);
			}
		}
		return result;
	}
	
	
	template< class MatrixType >
	const std::pair< MatrixType, MatrixType > Eig(const MatrixType &m)
	{
		//TODO with QR Decomposition
		return m;
	}
	

	
	
	
	template< class M, typename T > 
	const M MatrixFromFile(const std::string &path, const char colSeparator = '\t') 
	{
		std::ifstream file(path.c_str(), std::ifstream::in);
		
		
		static const int MAX_BUFF_SIZE = 1024;
		char buffer[MAX_BUFF_SIZE];
		
		if (!file.good()) {
			file.close();
			M m(1,1);
			
			std::cout << "File Not Found at: " << path << std::endl;
			return m;
		}
		
		int rows = 0, columns = 0;
		
		std::vector< std::vector< T > > tempMatrix;
		
		while (!file.eof()) {
			file.getline(buffer, MAX_BUFF_SIZE);
			const std::string line(buffer);
			tempMatrix.push_back(utils::SplitAndParseString< T >(line, colSeparator));
			if (tempMatrix[rows].empty()) {
				break;
			}
			
			columns = (int) tempMatrix[rows].size();
			
			NSE_SOLVER_ASSERT(columns == 0 || columns == (int) tempMatrix[rows].size());
			
			rows++;
		}
		
		
		M m(rows, columns);
		
		int i = 0;
		for (typename std::vector< std::vector< T > >::iterator row = tempMatrix.begin(); row != tempMatrix.end(); ++row) {
			int j = 0;
			
			for (typename std::vector< T >::iterator it = row->begin(); it != row->end(); ++it) {
				//std::cout << *it << ".\t";
				m(i,j++) = *it;
			}
			std::cout << std::endl;
			i++; 
		}
		
		
		file.close();
		return m;
	}
	
}

#endif

